#########################################################################
#  File-Name: pobe_replication.R
#  Date: June 27th 2024
#  Author: Robert Vidigal
#  Purpose: POBE Replication
#  Out File: pobe_replication.log
#  Data In: study1dynata.RData; study2fb.RData; study3bovitz.RData
#  Data Out: Analyses, Tables, and Figures.
#  Prev file: None
#  Status: On-going
#  Machine: Mac Book PRO (Silicon M1 Chip)
#########################################################################
rm(list=ls()); gc(); objects()
# To execute an R file from the command line (bash).
# % R CMD BATCH ~/pobe_replication.R ~/logs/pobe_replication.log

library(lavaan)
library(ltm)
library(KernSmoothIRT)
library(psychometric)
library(ggplot2);

# # -----------------------------------------------------------------------
# STUDY 1 (DYNATA) - SPLIT BALLOT EXPERIMENT
# # -----------------------------------------------------------------------
### 6-Items used.
### Foreign Aid it is the area that the U.S. Federal Government least spends money (T)
### At present, there are three female justices on the U.S. Supreme Court. (T)
### Natural gas is the leading source of electricity in the U.S. (T)
### The political office now held by Angela Merkel is President of the European Union. (F)
#### A majority vote (one-half) is required for the U.S. Senate and U.S. House to override a presidential veto. (F)
#### The Republicans have a majority of seats in the U.S. House of Representatives (F)
setwd("~/Library/Mobile Documents/com~apple~CloudDocs/PhD/ARTIGOS Trabalhando/")
load("POBE/data/study1dynata.RData")

# SUBSETS -----------------------------------------------------------------
pkMC1<-subset(df1, CondM2==1) #  # Current MC
pkMC1DK<-subset(df1, CondM2==2) # Current MC with DK
pkMC1both<-subset(df1, CondM2==1 | CondM2==2)
pkTF<-subset(df1, CondM2==5) # TF 
pkTFDK<-subset(df1, CondM2==6) # TF with DK
pkTFboth<-subset(df1, CondM2==5 | CondM2==6)
# Note: this experiment contained 2 more conditions involving
# the MC format with a pledge for not searching for answers.
# (not used here)

# ALPHAS  -------------------------------------------------------------------
pkalphaMC1<-subset(pkMC1, select=c("spending1", "energy1", "women1", 
                                  "merkel1", "veto1", "majority1"))
pkalphaMC1<<-na.omit(pkalphaMC1)

pkalphaTF1<-subset(pkTF, select=c("spending3", "energy3", "women3",  
                                 "merkel3", "veto3", "majority3"))
pkalphaTF1<-na.omit(pkalphaTF1)

alphapkMC1<-psych::alpha(pkalphaMC1, keys=NULL,cumulative=FALSE,title=NULL, max=10,
                      na.rm=TRUE,check.keys=TRUE,n.iter=1000,delete=TRUE); 
print(alphapkMC1) # .43

alphapkTF1<-psych::alpha(pkalphaTF1, keys=NULL,cumulative=FALSE,title=NULL, max=10,
                        na.rm=TRUE,check.keys=TRUE,n.iter=1000,delete=TRUE); 
print(alphapkTF1) # .69

# Dimensions
psych::fa.parallel(pkalphaMC1) # 1component
psych::fa.parallel(pkalphaTF1) # 2component

### UNCERTAINTY -----------------------------------------------------------
df1$uncertainty<-df1$spendingU+df1$womenU+df1$energyU+df1$merkelU+df1$majorityU+df1$vetoU
prop.table(table(df1$uncertainty>3)) # 72% 4 or more uncertain answers
psych::describe(df1$uncertainty)
#hist(df1$uncertainty) # 4.5 uncertain response on average

# # -----------------------------------------------------------------------
### FIGURE 1 - DK Differences Female - Male
# # -----------------------------------------------------------------------
# Subsets
male<-subset(df1, gender==1) # 1103
female<-subset(df1, gender==0) # 1219

# Gender Uncertainty
t.test(female$uncertainty, male$uncertainty) # Female less certain (.81, p<.001)

# Uncertainty Model
lmUNCERTAIN1<-lm(uncertainty~female+edu+female*edu, data=df1); summary(lmUNCERTAIN1)

malepkMC1DK<-subset(pkMC1DK, gender==1)
malepkTFDK<-subset(pkTFDK, gender==1)
femalepkMC1DK<-subset(pkMC1DK, gender==0)
femalepkTFDK<-subset(pkTFDK, gender==0)

tbl_gender_dk<-data.frame(
  "quest"=c(rep(c("Federal Spending", "U.S. Energy Source", "Women in SCOTUS", 
            "Presidential Veto", "Merkel Office", "U.S. House Majority"),2)),
  "format"=c(rep("Multiple-choice",6), rep("Belief certainty",6)),
  "dif"=c((psych::describe(femalepkMC1DK$spendingdk)[[3]] - psych::describe(malepkMC1DK$spendingdk)[[3]]),
          (psych::describe(femalepkMC1DK$energydk)[[3]] - psych::describe(malepkMC1DK$energydk)[[3]]),
          (psych::describe(femalepkMC1DK$womendk)[[3]] - psych::describe(malepkMC1DK$womendk)[[3]]),
          (psych::describe(femalepkMC1DK$vetodk)[[3]] - psych::describe(malepkMC1DK$vetodk)[[3]]),
          (psych::describe(femalepkMC1DK$merkeldk)[[3]] - psych::describe(malepkMC1DK$merkeldk)[[3]]),
          (psych::describe(femalepkMC1DK$majoritydk)[[3]] - psych::describe(malepkMC1DK$majoritydk)[[3]]),
          (psych::describe(femalepkTFDK$spendingdk)[[3]] - psych::describe(malepkTFDK$spendingdk)[[3]]),
          (psych::describe(femalepkTFDK$energydk)[[3]] - psych::describe(malepkTFDK$energydk)[[3]]),
          (psych::describe(femalepkTFDK$womendk)[[3]] - psych::describe(malepkTFDK$womendk)[[3]]),
          (psych::describe(femalepkTFDK$vetodk)[[3]] - psych::describe(malepkTFDK$vetodk)[[3]]),
          (psych::describe(femalepkTFDK$merkeldk)[[3]] - psych::describe(malepkTFDK$merkeldk)[[3]]),
          (psych::describe(femalepkTFDK$majoritydk)[[3]] - psych::describe(malepkTFDK$majoritydk)[[3]])),
  "lower"=c(t.test(femalepkMC1DK$spendingdk, malepkMC1DK$spendingdk)$conf.int[1],
            t.test(femalepkMC1DK$energydk, malepkMC1DK$energydk)$conf.int[1],
            t.test(femalepkMC1DK$womendk, malepkMC1DK$womendk)$conf.int[1],
            t.test(femalepkMC1DK$vetodk, malepkMC1DK$vetodk)$conf.int[1],
            t.test(femalepkMC1DK$merkeldk, malepkMC1DK$merkeldk)$conf.int[1],
            t.test(femalepkMC1DK$majoritydk, malepkMC1DK$majoritydk)$conf.int[1],
            t.test(femalepkTFDK$spendingdk, malepkTFDK$spendingdk)$conf.int[1],
            t.test(femalepkTFDK$energydk, malepkTFDK$energydk)$conf.int[1],
            t.test(femalepkTFDK$womendk, malepkTFDK$womendk)$conf.int[1],
            t.test(femalepkTFDK$vetodk, malepkTFDK$vetodk)$conf.int[1],
            t.test(femalepkTFDK$merkeldk, malepkTFDK$merkeldk)$conf.int[1],
            t.test(femalepkTFDK$majoritydk, malepkTFDK$majoritydk)$conf.int[1]),
  "upper"= c(t.test(femalepkMC1DK$spendingdk, malepkMC1DK$spendingdk)$conf.int[2],
             t.test(femalepkMC1DK$energydk, malepkMC1DK$energydk)$conf.int[2],
             t.test(femalepkMC1DK$womendk, malepkMC1DK$womendk)$conf.int[2],
             t.test(femalepkMC1DK$vetodk, malepkMC1DK$vetodk)$conf.int[2],
             t.test(femalepkMC1DK$merkeldk, malepkMC1DK$merkeldk)$conf.int[2],
             t.test(femalepkMC1DK$majoritydk, malepkMC1DK$majoritydk)$conf.int[2],
             t.test(femalepkTFDK$spendingdk, malepkTFDK$spendingdk)$conf.int[2],
             t.test(femalepkTFDK$energydk, malepkTFDK$energydk)$conf.int[2],
             t.test(femalepkTFDK$womendk, malepkTFDK$womendk)$conf.int[2],
             t.test(femalepkTFDK$vetodk, malepkTFDK$vetodk)$conf.int[2],
             t.test(femalepkTFDK$merkeldk, malepkTFDK$merkeldk)$conf.int[2],
             t.test(femalepkTFDK$majoritydk, malepkTFDK$majoritydk)$conf.int[2]))

tbl_gender_dk

plot_gender_dk<-ggplot(tbl_gender_dk, aes(x=reorder(quest, -dif), 
                                          y=dif, shape=format, color=format)) + 
  geom_errorbar(aes(ymin=lower, ymax=upper), 
                width=0, size=1, position=position_dodge(.2), na.rm = T) + 
  theme_classic() + geom_point(position=position_dodge(.2), na.rm = T, size=2) + 
  ylab("'Don't Know' Differences \n (Women – Men)") + xlab("") + 
  theme(legend.position="top", legend.title=element_blank(), 
        plot.title = element_text(size = 12, face = "bold"), 
        legend.text=element_text(size=12), 
        axis.text.x = element_text(size=10, face="bold"),  
        axis.title.y = element_text(size=12)) + 
  guides(col=guide_legend(byrow=T)) + 
  scale_x_discrete(labels=c("Angela Merkel's \n Office", "Women in \n SCOTUS", 
                            "U.S. Federal \n Spending (least)", "Presidential \n Veto override", 
                            "U.S. \n House \n Majority", "U.S. \n Energy \n Source"))  + 
  scale_y_continuous(limits=c(-.3,.3), breaks=c(-.3,-.2,-.1,0,.1,.2,.3)) + 
  scale_colour_manual( labels = c("Belief certainty", "Multiple-choice"), 
                       values=c("black", "grey66"))  +  
  geom_hline(yintercept=0, linetype="dashed")

jpeg("POBE/figures/Figure1.jpg", units="in", width = 6, height = 3, res=1200)
plot_gender_dk
dev.off()

t.test(femalepkMC1DK$pk_dk, femalepkTFDK$pk_dk)

# # --------------------------------------------------------------------------------------
# CHEATING: CATCH PAYNE (1875) & SELF-REPORTS --- NO DIFS
# # --------------------------------------------------------------------------------------
#prop.table(table(pkMC1$catch_payne==1875)) # Catch Q
#prop.table(table(pkMC1$cheat)) # Self-report
#
#prop.table(table(pkTF$catch_payne==1875)) # Catch Q
#prop.table(table(pkTF$cheat)) # Self-report

# -----------------------------------------------------------------------
# -----------------------------------------------------------------------
# -----------------------------------------------------------------------

# # -----------------------------------------------------------------------
# STUDY 2 (META ADS) - SPLIT-BALLOT EXPERIMENT
# # -----------------------------------------------------------------------
# know1 "Term Limits"
# know2 "Presidential Veto"
# know3 "U.S. House Majority"
# know4 "Women in SCOTUS"
# know5 "Midterm Elections Proportion"
# know6 "U.K. Prime Minister"

# -----------------------------------------------------------------------
# CSIP DATA - SPLIT BALLOT EXPERIMENT FOR SEARCHING BEHAVIOR
# # -----------------------------------------------------------------------
load("POBE/data/study2fb.RData")
df2<-df2 %>% mutate_at(vars(contains("know")), as.numeric)

### KNOW MC
# -----------------------------------------------------------------------
table(df2$know1MC)
table(df2$know2MC)
table(df2$know3MC)
table(df2$know4MC)
table(df2$know5MC)
table(df2$know6MC)
df2$knowMCscale<-df2$know1MC+df2$know2MC+df2$know3MC+df2$know4MC+df2$know5MC+df2$know6MC
table(df2$knowMCscale)

### SEARCH MC
# -----------------------------------------------------------------------
table(df2$know1MC_search)
psych::describe(df2$know1MC_timeSearch)

table(df2$know2MC_search)             
psych::describe(df2$know2MC_timeSearch)         

table(df2$know3MC_search)             
psych::describe(df2$know3MC_timeSearch)       

table(df2$know4MC_search)             
psych::describe(df2$know4MC_timeSearch)       

table(df2$know5MC_search)             
psych::describe(df2$know5MC_timeSearch)       

table(df2$know6MC_search)             
psych::describe(df2$know6MC_timeSearch)

table(df2$knowMCscale_search)

# UNCERTAINTY
# # -----------------------------------------------------------------------
table(df2$know1U)
table(df2$know2U)
table(df2$know3U)
table(df2$know4U)
table(df2$know5U)
table(df2$know6U)

# 57% gave at least 1 uncertain answer (6Q) 24% ZERO uncertain answers
prop.table(table(df2$uncertainty>1)) 

psych::describe(df2$uncertainty) # 2.04 mean

### KNOW TF
# # -----------------------------------------------------------------------
table(df2$know1TF)
table(df2$know2TF)
table(df2$know3TF)
table(df2$know4TF)
table(df2$know5TF)
table(df2$know6TF)
table(df2$knowTFscale)  

# UNCERTAINTY
# -----------------------------------------------------------------------
table(df2$uknow1TF)
table(df2$uknow2TF)
table(df2$uknow3TF)
table(df2$uknow4TF)
table(df2$uknow5TF)
table(df2$uknow6TF)

### SEARCH TF
# -----------------------------------------------------------------------
table(df2$know1TF_search)
psych::describe(df2$know1TF_timeSearch)

table(df2$know2TF_search)             
psych::describe(df2$know2TF_timeSearch)         

table(df2$know3TF_search)             
psych::describe(df2$know3TF_timeSearch)       

table(df2$know4TF_search)             
psych::describe(df2$know4TF_timeSearch)       

table(df2$know5TF_search)             
psych::describe(df2$know5TF_timeSearch)       

table(df2$know6TF_search)             
psych::describe(df2$know6TF_timeSearch)

table(df2$knowTFscale_search)

# # -----------------------------------------------------------------------
### DIFFERENCES BETWEEN MC and TF
# # -----------------------------------------------------------------------
# Means
t.test(df2$knowMCscale, df2$knowTFscale) # 2.61 vs 2.76
psych::describe(df2$knowMCscale); psych::describe(df2$knowTFscale)

# LESS INFORMATION SEARCH ON BELIEF CERTAINTY SCALE
t.test(df2$know1MC_search, df2$know1TF_search) # .001
t.test(df2$know2MC_search, df2$know2TF_search) # .003
t.test(df2$know3MC_search, df2$know3TF_search) # NO DIF .14
t.test(df2$know4MC_search, df2$know4TF_search) # .002
t.test(df2$know5MC_search, df2$know5TF_search) # .003
t.test(df2$know6MC_search, df2$know6TF_search) # .02

t.test(df2$knowMCscale_search, df2$knowTFscale_search) # .001, .35 vs 22

# # -----------------------------------------------------------------------
### FIGURE 2 - INFORMATION SEARCH 
# # -----------------------------------------------------------------------
tbl_search <- data.frame("quest"=c(rep(c("Term Limits", 
                                         "Presidential Veto", 
                                         "U.S. House Majority", 
                                         "Women in SCOTUS", 
                                         "Midterm Elections Proportion", 
                                         "U.K. Prime Minister"), 2)),
                  "mean"=c(psych::describe(df2$know1MC_search)[[3]], 
                           psych::describe(df2$know2MC_search)[[3]], 
                           psych::describe(df2$know3MC_search)[[3]], 
                           psych::describe(df2$know4MC_search)[[3]],
                           psych::describe(df2$know5MC_search)[[3]],
                           psych::describe(df2$know6MC_search)[[3]],
                           
                           psych::describe(df2$know1TF_search)[[3]], 
                           psych::describe(df2$know2TF_search)[[3]], 
                           psych::describe(df2$know3TF_search)[[3]], 
                           psych::describe(df2$know4TF_search)[[3]], 
                           psych::describe(df2$know5TF_search)[[3]],
                           psych::describe(df2$know6TF_search)[[3]]),
                  
                  "format"=c(rep("Multiple-choice",6), rep("Belief certainty",6)),
                  
                  "SE"=c(psych::describe(df2$know1MC_search)[[13]], 
                         psych::describe(df2$know2MC_search)[[13]], 
                         psych::describe(df2$know3MC_search)[[13]], 
                         psych::describe(df2$know4MC_search)[[13]],
                         psych::describe(df2$know5MC_search)[[13]],
                         psych::describe(df2$know6MC_search)[[13]],
                         
                         psych::describe(df2$know1TF_search)[[13]], 
                         psych::describe(df2$know2TF_search)[[13]], 
                         psych::describe(df2$know3TF_search)[[13]], 
                         psych::describe(df2$know4TF_search)[[13]], 
                         psych::describe(df2$know5TF_search)[[13]],
                         psych::describe(df2$know6TF_search)[[13]]))
tbl_search

tbl_search$quest<-factor(tbl_search$quest, levels=c("Term Limits", 
                                        "Presidential Veto", 
                                        "U.S. House Majority", 
                                        "Women in SCOTUS", 
                                        "Midterm Elections Proportion", 
                                        "U.K. Prime Minister"))
require(ggplot2)

jpeg("POBE/figures/Figure2.jpg", units="in", width = 10, height = 5, res=1200)

ggplot(tbl_search, aes(x=quest, y=mean, shape=format, colour=format)) + 
  geom_errorbar(aes(ymin=mean-1.64*SE, ymax=mean+1.64*SE), width=0,  size=1.5, 
                position=position_dodge(.2), na.rm = T) + 
  theme_classic() + geom_point(position=position_dodge(.2), na.rm = T, size=3) + 
  ylab("Answer look-up \n (screen-switching)") + 
  xlab("") + theme(legend.position="top", 
                   legend.title=element_blank(), 
                   legend.text=element_text(size=16), 
                   axis.text.x = element_text(face='bold', size=12), 
                   axis.text.y = element_text(face='bold', size=12), 
                   axis.title.y = element_text(size=16)) + 
  scale_x_discrete(labels=c("U.S. Congress \n Term Limits", 
                            "Presidential \n Veto", 
                            "U.S. House \n Majority", 
                            "Women in \n SCOTUS", 
                            "U.S. House/Senate \n seats in play", 
                            "U.K. \n Prime Minister")) + 
  scale_colour_manual(labels = c("Belief certainty", "Multiple-choice"), 
                      values=c("black", "grey66"))  + 
  scale_y_continuous(breaks=c(0,.02,.04,.06,.08,.1), labels = scales::percent) 

dev.off()

# # -----------------------------------------------------------------------
# KNOW SEARCH Scale Levels
# # -----------------------------------------------------------------------
df2$know1_search<-ifelse(is.na(df2$know1MC_search), df2$know1TF_search,df2$know1MC_search)
df2$know2_search<-ifelse(is.na(df2$know2MC_search), df2$know2TF_search,df2$know2MC_search)
df2$know3_search<-ifelse(is.na(df2$know3MC_search), df2$know3TF_search,df2$know3MC_search)
df2$know4_search<-ifelse(is.na(df2$know4MC_search), df2$know4TF_search,df2$know4MC_search)
df2$know5_search<-ifelse(is.na(df2$know5MC_search), df2$know5TF_search,df2$know5MC_search)
df2$know6_search<-ifelse(is.na(df2$know6MC_search), df2$know6TF_search,df2$know6MC_search)
df2$know_search<-df2$know1_search+df2$know2_search+df2$know3_search+
                 df2$know4_search+df2$know5_search+df2$know6_search

table(df2$know_search, df2$pkMC)

# means ttest
t.test(df2$knowMCscale_search, df2$knowTFscale_search) 
table(df2$knowMCscale_search>0); table(df2$knowTFscale_search>0); table(df2$pkMC) 

# Correlation screen-switch and knowledge score
summary(lm(knowMCscale~know_search, data=df2))
summary(lm(knowTFscale~know_search, data=df2))

cor.test(df2$knowMCscale, df2$know_search)
cor.test(df2$knowTFscale, df2$know_search)

# proportions chi-squared test
prop.test(x = c(144, 102), n = c(959, 961))

# ALPHAS  -------------------------------------------------------------------
pkalphaMC2<-subset(df2, select=c("know1MC", "know2MC", "know3MC", 
                                 "know4MC", "know5MC", "know6MC"))
pkalphaMC2<-na.omit(pkalphaMC2)

pkalphaTF2<-subset(df2, select=c("know1TF", "know2TF", "know3TF", 
                                 "know4TF", "know5TF", "know6TF"))
pkalphaTF2<-na.omit(pkalphaTF2)

alphapkMC2<-psych::alpha(pkalphaMC2, keys=NULL,cumulative=FALSE,title=NULL, max=10,
                        na.rm=TRUE,check.keys=TRUE,n.iter=1000,delete=TRUE); 
print(alphapkMC2) # .57

alphapkTF2<-psych::alpha(pkalphaTF2, keys=NULL,cumulative=FALSE,title=NULL, max=10,
                        na.rm=TRUE,check.keys=TRUE,n.iter=1000,delete=TRUE); 
print(alphapkTF2) # .71

# Dimensions
psych::fa.parallel(pkalphaMC2) # 1component
psych::fa.parallel(pkalphaTF2) # 2component

cor.test(df2$uncertainty, df2$knowTFscale)

# GENDER ANALYSIS
# # -----------------------------------------------------------------------
table(df2$gender)
df2$women<-ifelse(df2$gender==2,1,0)

male2<-subset(df2, gender==1)
female2<-subset(df2, gender==2)
t.test(male2$knowMCscale, female2$knowMCscale) # gap 3.25 vs 2.37
t.test(male2$knowTFscale, female2$knowTFscale) # gap 3.45 vs 2.36

t.test(male2$knowMCscale, male2$knowTFscale) 
t.test(female2$knowMCscale, female2$knowTFscale) 

# Gender Search
t.test(male2$know_search, female2$know_search) # NO DIF ON GENDER ON SEARCH p=.49

# Gender Uncertainty
t.test(male2$uncertainty, female2$uncertainty) # Women more uncertain (.47 p<.001)

# DiD
df2$knowscale<-ifelse(is.na(df2$knowMCscale), df2$knowTFscale, df2$knowMCscale)
table(is.na(df2$knowscale))
# pkMC==1 is multiple-choice
did2<-lm(knowscale~women+pkMC+pkMC*women, data=df2); summary(did2)
all.effects2 <- effects::allEffects(mod = did2)
summary(all.effects2)

# Uncertainty Model
lmUNCERTAIN<-lm(uncertainty~women+education1+women*education1, data=df2); summary(lmUNCERTAIN)

# # -----------------------------------------------------------------------
# SEARCHERS AND NON-SEARCHERS
# # -----------------------------------------------------------------------
searchers<-subset(df2, know_search!=0)
nonsearchers<-subset(df2, know_search==0)
t.test(searchers$uncertainty, nonsearchers$uncertainty) # p=.16

# # -----------------------------------------------------------------------
# DEVICE TYPE
# # -----------------------------------------------------------------------
device_desktop<-subset(df2, device_mobile==0)
device_phone<-subset(df2, device_mobile==1)
t.test(device_desktop$know_search, device_phone$know_search) # NO DIFFENCE IN SEARCHING p=.28

# # -----------------------------------------------------------------------
# BALANCE TESTS OF RANDOMIZATION FOR STUDY 2
# # -----------------------------------------------------------------------
df2MC<-subset(df2, pkMC==1); df2TF<-subset(df2, pkMC==0)
t.test(df2MC$women, df2TF$women) #ok
t.test(df2MC$education, df2TF$education) #ok
t.test(df2MC$income, df2TF$income) #ok

# # -----------------------------------------------------------------------
# # -----------------------------------------------------------------------
# # -----------------------------------------------------------------------

# -----------------------------------------------------------------------
# STUDY 3 (BOVITZ) - SPLIT BALLOT EXPERIMENT
# -----------------------------------------------------------------------
load("POBE/data/study3bovitz.RData")

### pk1 no term limits (TRUE)
### pk2 prime minister (FALSE)
### pk3 dems house majority (TRUE)
### pk4 3 women supreme court (FALSE)
### pk5 2/3 presidential veto (FALSE)

# PK TREATMENT VARIABLE
table(df3$pk) # 0 MC and 1 TF

pkMCdf3<-subset(df3, pk==0)
pkTFdf3<-subset(df3, pk==1)

# No difference in self-report cheating
t.test(pkMCdf3$cheating, pkTFdf3$cheating) # p=.79

# # -----------------------------------------------------------------------
## ALPHA
# # -----------------------------------------------------------------------
pkalphaMC3<-subset(pkMCdf3, select=c("pk1MC", "pk2MC", "pk3MC", "pk4MC", "pk5MC"))
pkalphaTF3<-subset(pkTFdf3, select=c("pk1TF", "pk2TF", "pk3TF",  "pk4TF", "pk5TF"))

pkalphaMC3<-na.omit(pkalphaMC3)
pkalphaTF3<-na.omit(pkalphaTF3)
alphapkMC3<-psych::alpha(pkalphaMC3, keys=NULL,cumulative=FALSE,title=NULL, 
                         max=10, na.rm=TRUE,check.keys=TRUE,n.iter=1000,delete=TRUE); 
alphapkMC3 # alpha .37
alphapkTF3<-psych::alpha(pkalphaTF3, keys=NULL,cumulative=FALSE,title=NULL, 
                         max=10, na.rm=TRUE,check.keys=TRUE,n.iter=1000,delete=TRUE); 
alphapkTF3 # alpha .66

# SCALES' ANALYSIS
# # ----
psych::describe(df3$pkMC); psych::describe(df3$pkTF)
hist(df3$pkMC)
hist(df3$pkTF)
dev.off()
t.test(df3$pkMC, df3$pkTF)

table(df3$know) # One Scale for both
table(is.na(df3$know)) # 66 NAs

# Dimensions
psych::fa.parallel(pkalphaMC3) # 1component
psych::fa.parallel(pkalphaTF3) # # 1 component

# UNCERTAINTY -------------------------------------------------------------
prop.table(table(df3$pkTFu>1)) 
# 44% gave 3+ uncertain answers (out of 5Q), 21% 0 Uncertain Answers

psych::describe(df3$pkTFu)
# 2.22 uncertain response on average

# # -----------------------------------------------------------------------
# DATA ANALYSIS
# # -----------------------------------------------------------------------

# OUTCOMES
# # ----

### TOLERANCE (Higher Values indicate Higher Tolerance)
# # ----
table(df3$tolerance)
pkalphaTOL<-subset(df3, select=c("tol1", "tol2", "tol3", "tol4", "tol5"))
pkalphaTOL<-na.omit(pkalphaTOL)
alphapkTOL<-psych::alpha(pkalphaTOL, keys=NULL,cumulative=FALSE,title=NULL, 
                         max=10, na.rm=TRUE,check.keys=TRUE,n.iter=1000,delete=TRUE); 
alphapkTOL # alpha .51

cor(df3$pkTFb, df3$tolerance, use="complete.obs") # Dichotomous
cor(df3$pkTF, df3$tolerance, use="complete.obs") # Polytomous stronger cor

### ENGAGEMENT (Higher Values indicate Higher Participation)
# # ----
table(df3$participation)
pkalphaPART<-subset(df3, select=c("engage_1", "engage_2", "engage_3", "engage_4", "engage_5"))
pkalphaPART<-na.omit(pkalphaPART)
alphapkPART<-psych::alpha(pkalphaPART, keys=NULL,cumulative=FALSE,title=NULL, 
                          max=10, na.rm=TRUE,check.keys=TRUE,n.iter=1000,delete=TRUE); 
alphapkPART # alpha .62

cor(df3$pkTFb, df3$participation, use="complete.obs") # Dichotomous
cor(df3$pkTF, df3$participation, use="complete.obs") # Polytomous stronger cor

### VACCINATION
# # ----
table(df3$vaccstatus)    
table(df3$vaccstatus2) # Higher values indicate endorsement of 2nd dose.

### Vaccine Beliefs (Higher values indicate pro vaccination views)
table(df3$vaccbeliefs)

pkalphaVACC<-subset(df3, select=c("vaccsafety", "vacceffective", "vaccencourage"))
pkalphaVACC<-na.omit(pkalphaVACC)
alphapkVACC<-psych::alpha(pkalphaVACC, keys=NULL,cumulative=FALSE,title=NULL, 
                          max=10, na.rm=TRUE,check.keys=TRUE,n.iter=1000,delete=TRUE); 
alphapkVACC # alpha .85

cor(df3$pkTFb, df3$vaccbeliefs, use="complete.obs") # Dichotomous
cor(df3$pkTF, df3$vaccbeliefs, use="complete.obs") # Polytomous stronger cor

### Vaccine Misinformation (Higher values indicate more misinformation beliefs)
table(df3$vaccmisinfo)

pkalphaMISINFO<-subset(df3, select=c("vaccmisinfo_1", "vaccmisinfo_2", "vaccmisinfo_3"))
pkalphaMISINFO<-na.omit(pkalphaMISINFO)
alphapkMISINFO<-psych::alpha(pkalphaMISINFO, keys=NULL,cumulative=FALSE,title=NULL, 
                             max=10, na.rm=TRUE,check.keys=TRUE,n.iter=1000,delete=TRUE); 
alphapkMISINFO # alpha .85

cor(df3$pkTFb, df3$vaccmisinfo, use="complete.obs") # Dichotomous
cor(df3$pkTF, df3$vaccmisinfo, use="complete.obs") # Polytomous stronger cor

# Controls & Others
# # ----
# Income
table(df3$income)

# Cheating
table(df3$cheating)

## PID (1 Rep, 2 Dem, 3 Ind)
table(df3$pid)
table(df3$demstrength)
table(df3$repstrength)
table(df3$ind_lean)

df3$partyid<-factor(df3$pid, levels=c(3, 1, 2), 
                    labels=c("Independent", "Republican", "Democrat"))
table(df3$partyid) # Independent as category base

df3$partyid2<-factor(df3$pid, levels=c(1, 2), labels=c("Republican", "Democrat"))
table(df3$partyid2); df3$partyid2 <- sjlabelled::set_label(df3$partyid2, label = "Party ID")

table(df3$partyid3) # Democrat as the category base

df3$pidstr<-factor(df3$pidst, levels=c(1,7), labels=c("Republicans", "Democrats"))
table(df3$pidstr)
df3$pidstr <- sjlabelled::set_label(df3$pidstr, label = "Party ID")

## Education (higher values higher education)
table(df3$education)

## Age
table(df3$age)

## Income
table(df3$income)

# # -----------------------------------------------------------------------
# MARGINAL EFFECTS MODELS
# # -----------------------------------------------------------------------
## TOLERANCE OUTCOME
reg_tolMC<-lm(tolerance~pkMC+education+female+age+income+partyid2+partyid2*pkMC, 
              data=df3); summary(reg_tolMC)
reg_tolTF<-lm(tolerance~pkTF+education+female+age+income+partyid2+partyid2*pkTF, 
              data=df3); summary(reg_tolTF)

margtolMC<-summary(margins::margins(reg_tolMC, terms=c("pkMC")))
margtolTF<-summary(margins::margins(reg_tolTF, terms=c("pkTF")))

## ENGAGEMENT OUTCOME
reg_partMC<-lm(participation~pkMC+education+female+age+income+partyid2+partyid2*pkMC, 
               data=df3); summary(reg_partMC)
reg_partTF<-lm(participation~pkTF+education+female+age+income+partyid2+partyid2*pkTF, 
               data=df3); summary(reg_partTF)

margpartMC<-summary(margins::margins(reg_partMC, terms=c("pkMC")))
margpartTF<-summary(margins::margins(reg_partTF, terms=c("pkTF")))

## VACCINATION BELIEFS
reg_VACCbeliefsMC<-lm(vaccbeliefs~pkMC+education+female+age+income+partyid2+partyid2*pkMC, 
                      data=df3); summary(reg_VACCbeliefsMC)
reg_VACCbeliefsTF<-lm(vaccbeliefs~pkTF+education+female+age+income+partyid2+partyid2*pkTF, 
                      data=df3); summary(reg_VACCbeliefsTF)

margVACCbeliefsMC<-summary(margins::margins(reg_VACCbeliefsMC, terms=c("pkMC")))
margVACCbeliefsTF<-summary(margins::margins(reg_VACCbeliefsTF, terms=c("pkTF")))


## VACCINE MISINFORMATION
reg_VACCmisinfoMC<-lm(vaccmisinfo~pkMC+education+female+age+income+partyid2+partyid2*pkMC, 
                      data=df3); summary(reg_VACCmisinfoMC)
reg_VACCmisinfoTF<-lm(vaccmisinfo~pkTF+education+female+age+income+partyid2+partyid2*pkTF, 
                      data=df3); summary(reg_VACCmisinfoTF)

margVACCmisinfoMC<-summary(margins::margins(reg_VACCmisinfoMC, terms=c("pkMC")))
margVACCmisinfoTF<-summary(margins::margins(reg_VACCmisinfoTF, terms=c("pkTF")))

# # -----------------------------------------------------------------------
### FIGURE 3 - PLOT AVG MARGINAL EFFECTS ACROSS PK TYPES
# # -----------------------------------------------------------------------
tblmarg<-data.frame(
  "cond"=c(rep("Political Tolerance",2),rep("Political Engagement",2),
           rep("Vaccine Beliefs",2),rep("Vaccine Misinfo",2)),
                    "knowtype"=c(rep(c("Multiple-choice", "Belief certainty"),4)), 
                    "ame"=c(margtolMC[6,2], margtolTF[6,2],
                            margpartMC[6,2], margpartTF[6,2], 
                            margVACCbeliefsMC[6,2], margVACCbeliefsTF[6,2], 
                            margVACCmisinfoMC[6,2], margVACCmisinfoTF[6,2]), 
                    "SE"=c(margtolMC[6,3], margtolTF[6,3],
                           margpartMC[6,3], margpartTF[6,3], 
                           margVACCbeliefsMC[6,3], margVACCbeliefsTF[6,3], 
                           margVACCmisinfoMC[6,3], margVACCmisinfoTF[6,3]))

tblmarg

tblmarg$cond<-factor(tblmarg$cond, levels=c("Political Tolerance", 
                                            "Political Engagement", 
                                            "Vaccine Beliefs", 
                                            "Vaccine Misinfo"), ordered=T)

require(ggplot2)
plotmarg2<-ggplot(tblmarg, aes(x=cond, y=ame, shape=knowtype, colour=knowtype)) + 
  geom_errorbar(aes(ymin=ame-1.96*SE, ymax=ame+1.96*SE), width=0, size=1, 
                position=position_dodge(.3), na.rm = T) +
  theme_classic() + geom_point(position=position_dodge(.3), na.rm = T, size=3) +
  theme(legend.position="top", legend.title=element_blank(), 
        plot.title = element_text(size = 12, face = "bold"), 
        legend.text=element_text(size=10)) + 
  geom_hline(yintercept = 0, linetype="dashed", alpha=.3) + 
  theme(axis.text.x = element_text(face='bold', size=10), 
        axis.title.y = element_text(size=12), 
        legend.text = element_text(size=14), 
        axis.text.y = element_text(face='bold', size=12)) + 
  xlab("") + 
  scale_x_discrete(labels=c("Political \n Tolerance", 
                            "Political \n Engagement", 
                            "Vaccine \n Beliefs", 
                            "Vaccine \n Misinformation")) + 
  ylim(-.1,.1) + ylab("Coefficient Estimate") +  guides(col=guide_legend(byrow=T)) + 
  scale_colour_manual(labels = c("Belief certainty", "Multiple-choice"), 
                      values=c("black", "grey66"))

jpeg("POBE/figures/Figure3.jpg", units="in", width = 6, height = 4, res=1200)

plotmarg2 

dev.off()

# # -----------------------------------------------------------------------
# FIGURE 4 - CORRECT-RESPONSE RATE RELIABILITY
# # -----------------------------------------------------------------------
# Study 1 WomenSC HouseMAJOR PresVETO
# Study 2 know1=TermLIMITS know2=PresVETO know3=HouseMAJOR know4=WomenSC
# Study 3 pk1=TermLIMITS pk3=HouseMAJOR pk4=WomenSC pk5=PresVETO

# The largest differences between studies are taken, unless they differ in question wording.

# WOMEN SUPREME COURT
# # -----------------------------------------------------------------------

# MC
prop.table(table(df1$women1)) # 43%
prop.table(table(df2$know4MC)) # 12%
prop.table(table(df3$pk4MC)) # 34%
# 30% variation
round(abs(prop.table(table(df1$women1))[[2]] - prop.table(table(df2$know4MC))[[2]]), 2)*100

# BC
prop.table(table(df1$women3)) # 63%
prop.table(table(df2$know4TF)) # 58% (15+43)
prop.table(table(df3$pk4TF)) # 52% (19+33)
# 11% variation (dif T/F orientation) or 6% variation (same wording)
round(abs(prop.table(table(df2$know4TF!=0))[[2]] - prop.table(table(df3$pk4TF!=0))[[2]]), 2)*100

# HOUSE MAJORITY
# # -----------------------------------------------------------------------
# MC
prop.table(table(df1$majority1)) # 58%
prop.table(table(df2$know3MC)) # 71%
prop.table(table(df3$pk3MC)) # 77%
# 19% variation
round(abs(prop.table(table(df1$majority1))[[2]] - prop.table(table(df3$pk3MC))[[2]]), 2)*100

# BC
prop.table(table(df1$majority3)) # 49% (10 + 39)
prop.table(table(df2$know3TF)) # 76% (26+50)
prop.table(table(df3$pk3TF)) # 84% (30+54)
# 35% variation (dif T/F) or 8% (same wording)
round(abs(prop.table(table(df2$know3TF!=0))[[2]] - prop.table(table(df3$pk3TF!=0))[[2]]), 2)*100

# TERM LIMITS
# # -----------------------------------------------------------------------
# MC
prop.table(table(df2$know1MC)) # 38%
prop.table(table(df3$pk1MC)) # 52%
# 14% variation
round(abs(prop.table(table(df2$know1MC))[[2]] -  prop.table(table(df3$pk1MC))), 2)*100

# BC
prop.table(table(df2$know1TF!=0)) # 57% (22+35)
prop.table(table(df3$pk1TF!=0)) # 56% (20+35)
# 1% variation
round(abs(prop.table(table(df2$know1TF!=0))[[2]] - prop.table(table(df3$pk1TF!=0))[[2]]), 2)*100

# UK PRIME MINISTER
# # -----------------------------------------------------------------------
# MC
prop.table(table(df2$know6MC)) # 56%
prop.table(table(df3$pk2MC)) # 71%
# 15% variation
round(abs(prop.table(table(df2$know6MC))[[2]] - prop.table(table(df3$pk2MC))[[2]]), 2)*100

# BC
prop.table(table(df2$know6TF)) # 60% (9+51)
prop.table(table(df3$pk2TF)) # 69% (13+56)
# 9% variation
round(abs(prop.table(table(df2$know6TF!=0))[[2]] - prop.table(table(df3$pk2TF!=0))[[2]]), 2)*100

# PRESIDENTIAL VETO
# # -----------------------------------------------------------------------
# MC
prop.table(table(df1$veto1)) # 61%
prop.table(table(df2$know2MC)) # 66%
prop.table(table(df3$pk5MC)) # 77%
# 16% variation
round(abs(prop.table(table(df1$veto1))[[2]] - prop.table(table(df3$pk5MC))[[2]]), 2)*100

# BC
prop.table(table(df1$veto3)) # 43% (15 + 28)
prop.table(table(df2$know2TF)) # 37% (13+24)
prop.table(table(df3$pk5TF)) # 46% (18+28)
# 10% variation
round(abs(prop.table(table(df2$know2TF!=0))[[2]] - prop.table(table(df3$pk5TF!=0))[[2]]), 2)*100

# # -----------------------------------------------------------------------
# PLOTTING RELIABILITY CORRECT-RESPONSE RATE
# # -----------------------------------------------------------------------
tblcr <- data.frame(
  "quest"=c(rep(c("U.S. Term Limits", 
                  "U.K. Prime Minister",
                  "Presidential Veto", 
                  "Women SCOTUS", 
                  "U.S. House Majority"), 2)),
            study=c(rep("Variation (∆)", 10)),
  "format"=c(rep("Multiple-choice", 5), rep("Belief certainty", 5)),
  "mean"=c(
    #MC
    round(abs(prop.table(table(df2$know1MC))[[2]] -  
                prop.table(table(df3$pk1MC)))[[2]], 2),
    round(abs(prop.table(table(df2$know6MC))[[2]] - 
                prop.table(table(df3$pk2MC))[[2]]), 2),
    round(abs(prop.table(table(df1$veto1))[[2]] - 
                prop.table(table(df3$pk5MC))[[2]]), 2),
    round(abs(prop.table(table(df1$women1))[[2]] - 
                prop.table(table(df2$know4MC))[[2]]), 2),
    round(abs(prop.table(table(df1$majority1))[[2]] - 
                prop.table(table(df3$pk3MC))[[2]]), 2),
    #BC
    round(abs(prop.table(table(df2$know1TF!=0))[[2]] - 
                prop.table(table(df3$pk1TF!=0))[[2]]), 2), 
    round(abs(prop.table(table(df2$know6TF!=0))[[2]] - 
                prop.table(table(df3$pk2TF!=0))[[2]]), 2),
    round(abs(prop.table(table(df2$know2TF!=0))[[2]] - 
                prop.table(table(df3$pk5TF!=0))[[2]]), 2),
    round(abs(prop.table(table(df2$know4TF!=0))[[2]] - 
                prop.table(table(df3$pk4TF!=0))[[2]]), 2),
    round(abs(prop.table(table(df2$know3TF!=0))[[2]] - 
                prop.table(table(df3$pk3TF!=0))[[2]]), 2)))

tblcr$format<-factor(tblcr$format, levels=c("Belief certainty", "Multiple-choice"))

tblcr$quest<-factor(tblcr$quest, levels=c("U.S. Term Limits", 
                                          "U.K. Prime Minister",
                                          "Presidential Veto", 
                                          "Women SCOTUS", 
                                          "U.S. House Majority"))
tblcr

# PLOTTING
# # --------------------------------------------------------------------------------------
jpeg("POBE/figures/Figure4.jpg", units="in", width = 10, height = 7, res=1200)
ggplot(tblcr, aes(x=quest, y=mean, fill=format)) + 
  theme_classic() + geom_bar(stat = "identity", width=.6,  
                             position =position_dodge(width = .65),
                             lwd=0.5, color="black") +
  ylab("Variation (∆) in correct-response rate") +
  xlab("") + theme(legend.position="top", 
                   legend.title=element_blank(), 
                   legend.text=element_text(size=16), 
                   axis.text.x = element_text(face='bold', size=12), 
                   axis.text.y = element_text(face='bold', size=12), 
                   axis.title.y = element_text(size=14)) +
  geom_text(aes(label=paste0(round(mean*100), "%")),
            color="white",
            size=4, 
            fontface = "bold",
            vjust = 1.5, 
            position = position_dodge(.65)) +
  scale_x_discrete(labels=c("U.S. Term Limits \n (S2, S3)", 
                            "U.K. Prime Minister \n (S2, S3)", 
                            "Presidential Veto \n (S1, S2, S3)", 
                            "Women in SCOTUS \n (S1, S2, S3)", 
                            "U.S. House Majority \n (S1, S2, S3)")) + 
  scale_fill_manual(labels = c("Belief certainty", "Multiple-choice"), 
                    values=c("black", "grey66"))  + 
  scale_y_continuous(breaks=c(0,.10,.20,.30,.40,.50), limits=c(0,.31), 
                     labels = scales::percent, expand = c(0, 0)) 
dev.off()

# # -----------------------------------------------------------------------
# # -----------------------------------------------------------------------
# # -----------------------------------------------------------------------

# # -----------------------------------------------------------------------
# MEASUREMENT MODELS (STUDY 1) - WITH DK 
# # -----------------------------------------------------------------------
MCss<-subset(df1, CondM2==1 | CondM2==2, #same results without DK
               select=c("spending1", "energy1", "women1", 
                        "veto1", "merkel1", "majority1", "female"))
MCss1<-subset(MCss, select=-c(female))

TFss<-subset(df1, CondM2==5 | CondM2==6, #same results without DK
              select=c("spending3", "energy3", "women3", 
                       "veto3", "merkel3", "majority3", "female"))
TFss1<-subset(TFss, select=-c(female))

# Recoding (scoring does not change results, but necessary for mirt)
TFss$spending3<-car::recode(TFss$spending3, "1 = 2; 0.5 = 1"); table(TFss$spending3)
TFss$energy3<-car::recode(TFss$energy3, "1 = 2; 0.5 = 1"); table(TFss$energy3)
TFss$women3<-car::recode(TFss$women3, "1 = 2; 0.5 = 1"); table(TFss$women3)
TFss$merkel3<-car::recode(TFss$merkel3, "1 = 2; 0.5 = 1"); table(TFss$merkel3)
TFss$veto3<-car::recode(TFss$veto3, "1 = 2; 0.5 = 1"); table(TFss$veto3)
TFss$majority3<-car::recode(TFss$majority3, "1 = 2; 0.5 = 1"); table(TFss$majority3)

# Renaming for plotting
TFss <- TFss %>% rename("Foreign Aid Spending"=spending3,
                      "U.S. Source of Energy"=energy3,
                      "Women in SCOTUS"=women3,
                      "Presidential Veto"=veto3,
                      "Angela Merkel's Office"=merkel3,
                      "U.S. House Majority"=majority3)

# Nominal Measurement Model
require(mirt); mirt.model1<-'polknow=1-6'

# --- MC
nomMC1<-mirt(MCss1, mirt.model1, itemtype="nominal", SE=T); 
round(coef(nomMC1, IRTpars=F,simplify=T)$items[1:6,c(1,5)], 3)
# The values is the first column (a1) reflect the item slopes, 
# while the values is the second column (d1) correspond to the item intercept. 

# OCC
plot(nomMC1, type = 'trace', which.items = c(1:6),  theta_lim = c(-3, 3), 
     main = "\n", par.settings = simpleTheme(lty=c(1,2,1),lwd=3),
                   points=FALSE,lines=TRUE, columns=3)

# --- BC
nomTF1<-mirt(TFss, mirt.model1, itemtype="nominal", SE=T); 
round(coef(nomTF1, IRTpars=F,simplify=T)$items[1:6,], 3)

# OCC
jpeg("POBE/figures/Appendix_OCC1_scoring.jpg", units="in", width = 8, height = 6, res=1200)
plot(nomTF1, type = 'trace', which.items = c(1:6),  theta_lim = c(-3, 3), 
     main = "\n", par.settings = simpleTheme(lty=c(2,1,1),
                                             lwd=3, 
                                             col.line = c("skyblue", "orange", "darkgreen")),
     auto.key=list(text=c("Wrong Answer/DK", "Uncertain Answer", "Certain Answer"),
                   points=FALSE,lines=TRUE, columns=3))
dev.off()

# --- WITH DK (SAME RESULTS CUZ DK SCORED AS 0)
# # --------------------------------------------------------------------------------------
#MCssDK<-subset(df1, CondM2==2, 
#               select=c("spending1", "energy1", "women1", 
#                        "veto1", "merkel1", "majority1"))
#TFssDK<-subset(df1, CondM2==6, 
#              select=c("spending3", "energy3", "women3", 
#                       "veto3", "merkel3", "majority3"))
## Recoding (scoring does not change results, but necessary for mirt)
#TFssDK$spending3<-car::recode(TFssDK$spending3, "1 = 2; 0.5 = 1"); table(TFssDK$spending3)
#TFssDK$energy3<-car::recode(TFssDK$energy3, "1 = 2; 0.5 = 1"); table(TFssDK$energy3)
#TFssDK$women3<-car::recode(TFssDK$women3, "1 = 2; 0.5 = 1"); table(TFssDK$women3)
#TFssDK$merkel3<-car::recode(TFssDK$merkel3, "1 = 2; 0.5 = 1"); table(TFssDK$merkel3)
#TFssDK$veto3<-car::recode(TFssDK$veto3, "1 = 2; 0.5 = 1"); table(TFssDK$veto3)
#TFssDK$majority3<-car::recode(TFssDK$majority3, "1 = 2; 0.5 = 1"); table(TFssDK$majority3)
#
#nomMC1DK<-mirt(MCssDK, mirt.model1, itemtype="nominal", SE=T); 
#nomTF1DK<-mirt(TFssDK, mirt.model1, itemtype="nominal", SE=T); 
#
#plot(nomTF1DK, type = 'trace', which.items = c(1:6),  theta_lim = c(-3, 3), 
#     main = "\n", par.settings = simpleTheme(lty=c(1,2,1),
#                                             lwd=3, 
#                                             col.line = c("skyblue", "orange", "darkgreen")),
#     auto.key=list(text=c("Uncertain Answer","Wrong Answer", "Certain Answer"),
#                   points=FALSE,lines=TRUE, columns=3))

# # -----------------------------------------------------------------------
# MEASUREMENT MODELS (STUDY 2)
# # -----------------------------------------------------------------------
MC2ss<-subset(df2, select=c("know1MC", "know2MC", "know3MC", 
                                  "know4MC", "know5MC", "know6MC", "women"))
MC2ss<-na.omit(MC2ss)

TF2ss<-subset(df2, select=c("know1TF", "know2TF", "know3TF", 
                             "know4TF", "know5TF", "know6TF", "women"))
TF2ss<-na.omit(TF2ss)

# Recoding (scoring does not change results, but necessary for mirt)
TF2ss$know1TF<-car::recode(TF2ss$know1TF, "1 = 2; 0.5 = 1"); table(TF2ss$know1MC)
TF2ss$know2TF<-car::recode(TF2ss$know2TF, "1 = 2; 0.5 = 1"); table(TF2ss$know2MC)
TF2ss$know3TF<-car::recode(TF2ss$know3TF, "1 = 2; 0.5 = 1"); table(TF2ss$know3MC)
TF2ss$know4TF<-car::recode(TF2ss$know4TF, "1 = 2; 0.5 = 1"); table(TF2ss$know4MC)
TF2ss$know5TF<-car::recode(TF2ss$know5TF, "1 = 2; 0.5 = 1"); table(TF2ss$know5MC)
TF2ss$know6TF<-car::recode(TF2ss$know6TF, "1 = 2; 0.5 = 1"); table(TF2ss$know6MC)

TF2ss <- TF2ss %>% rename("U.S. Terms Limits"=know1TF,
                      "Presidential Veto"=know2TF,
                      "U.S. House Majority"=know3TF,
                      "Women in SCOTUS"=know4TF,
                      "2022 Election Seats"=know5TF,
                      "U.K. Prime Minister"=know6TF)

# Nominal
require(mirt); mirt.model2<-'know=1-6'

# MC
nomMC2<-mirt(MC2ss, mirt.model2, itemtype="nominal", SE=T); 
round(coef(nomMC2, simplify=T)$items[1:6,c(1,5)], 3)

# OCC
plot(nomMC2, type = 'trace', which.items = c(1:6),  theta_lim = c(-3, 3), 
     main = "\n", par.settings = simpleTheme(lty=c(1),lwd=3),
     auto.key=list(text=c(""),
                   points=FALSE,lines=TRUE, columns=3))

# BC
nomTF2<-mirt(TF2ss, mirt.model2, itemtype="nominal", SE=T);
round(coef(nomTF2, simplify=T)$items[1:6,], 3)

# OCC
jpeg("POBE/figures/Appendix_OCC2_scoring.jpg", units="in", width = 8, height = 6, res=1200)
plot(nomTF2, type = 'trace', which.items = c(1:6),  theta_lim = c(-3, 3), 
     main = "\n", par.settings = simpleTheme(lty=c(2,1,1),
                                             lwd=3, 
                                             col.line = c("skyblue", "orange", "darkgreen")),
     auto.key=list(text=c("Wrong Answer", "Uncertain Answer", "Certain Answer"),
                   points=FALSE,lines=TRUE, columns=3))
dev.off()

# # -----------------------------------------------------------------------
## MEASUREMENT MODELS (STUDY 3)
# # -----------------------------------------------------------------------
MC3ss<-subset(pkMCdf3, select=c("pk1MC", "pk2MC", "pk3MC", "pk4MC", "pk5MC", "female"))
MC3ss<-na.omit(MC3ss)

TF3ss<-subset(pkTFdf3, select=c("pk1TF", "pk2TF", "pk3TF",  "pk4TF", "pk5TF", "female"))
TF3ss<-na.omit(TF3ss)

# Recoding (scoring does not change results, but necessary for mirt)
TF3ss$pk1TF<-car::recode(TF3ss$pk1TF, "1 = 2; 0.5 = 1"); table(TF3ss$pk1TF)
TF3ss$pk2TF<-car::recode(TF3ss$pk2TF, "1 = 2; 0.5 = 1"); table(TF3ss$pk2TF)
TF3ss$pk3TF<-car::recode(TF3ss$pk3TF, "1 = 2; 0.5 = 1"); table(TF3ss$pk3TF)
TF3ss$pk4TF<-car::recode(TF3ss$pk4TF, "1 = 2; 0.5 = 1"); table(TF3ss$pk4TF)
TF3ss$pk5TF<-car::recode(TF3ss$pk5TF, "1 = 2; 0.5 = 1"); table(TF3ss$pk5TF)

TF3ss <- TF3ss %>% rename("Terms Limits House"=pk1TF,
                      "U.K. Prime Minister"=pk2TF,
                      "U.S. House Majority"=pk3TF,
                      "Women in SCOTUS"=pk4TF,
                      "Presidential Veto"=pk5TF)

# mirt
mirt.model3<-'polknow=1-5'

# MC
nomMC3<-mirt(MC3ss, mirt.model3, itemtype="nominal", SE=T); nomMC3; 
round(coef(nomMC3, IRTpars=F, simplify=T)$items[1:5,c(1,5)], 3)

plot(nomMC3, type = 'trace', which.items = c(1:5),  theta_lim = c(-3, 3), 
     main = "\n", par.settings = simpleTheme(lty=c(1,2,1),
                                             lwd=3))
# TF
nomTF3<-mirt(TF3ss, mirt.model3, itemtype="nominal", SE=T); nomTF3; 
round(coef(nomTF3, IRTpars=F, simplify=T)$items[1:5,], 3)

jpeg("POBE/figures/Appendix_OCC3_scoring.jpg", units="in", width = 8, height = 6, res=1200)
plot(nomTF3, type = 'trace', which.items = c(1:5),  theta_lim = c(-3, 3), 
     main = "\n", par.settings = simpleTheme(lty=c(2,1,1),
                                             lwd=3, 
                                             col.line = c("skyblue", "orange", "darkgreen")),
     auto.key=list(text=c("Wrong Answer","Uncertain Answer", "Certain Answer"),
                   points=FALSE,lines=TRUE, columns=3))
dev.off()
# # -----------------------------------------------------------------------
# FIGURE 5  - INFORMATION ANALYSIS (TIF)
# # -----------------------------------------------------------------------
set.seed(10012)
jpeg("POBE/figures/Figure5.jpg", units="in", height = 4, width = 12, res=1200)
layout(matrix(1:3, ncol = 3, byrow = TRUE))

# # -----------------------------------------------------------------------
### TIF STUDY 1
# # -----------------------------------------------------------------------
# The TIF is the sum of the item information functions (IIF) in a test, 
# and indicates the precision of measurement that can be achieved with 
# the test at any value of the latent variable.

# Test Information 
Theta1 <- matrix(seq(-3,3,.05))
plot(Theta1, testinfo(nomMC1, Theta1), type = 'l', main = 'Study 1 (Dynata)', 
     ylab = "", xlab=substitute(paste('θ')), ylim=c(0,12), lwd=3, col = 'gray45', lty=1,
     cex.main=2, cex.axis=2, cex.lab=3, frame.plot=F, mgp = c(3.5, 0.5, -0.5))
lines(Theta1, testinfo(nomTF1, Theta1), col = 'black', lty=1, lwd=3)
mtext(side=2, line=2.5, "Information", cex=1.5, font=1)
legend("topright", legend=c("Belief certainty", "Multiple-choice"),
       col=c("black", "gray45"), lty=c(1,1), lwd=4, cex=2, bty="n")

# Area under the curves.
areainfo(nomMC1, c(-3,3)) # 5.28
areainfo(nomTF1, c(-3,3)) # 12.53

# # -----------------------------------------------------------------------
### TIF STUDY 2 
# # -----------------------------------------------------------------------
# The TIF is the sum of the item information functions (IIF) in a test, 
# and indicates the precision of measurement that can be achieved with the 
# test at any value of the latent variable
Theta2 <- matrix(seq(-3,3,.05))
plot(Theta2, testinfo(nomMC2, Theta2), type = 'l', main = 'Study 2 (Meta Ads)', 
     ylab = "", xlab=substitute(paste('θ')), ylim=c(0,12), lwd=3, col = 'gray40', lty=1, 
     cex.main=2, cex.axis=2, cex.lab=3, frame.plot=F, mgp = c(3.5, 0.5, -0.5))
lines(Theta2, testinfo(nomTF2, Theta2), col = 'black', lty=1, lwd=3)
mtext(side=2, line=2.5, "Information", cex=1.5, font=1)
legend("topright", legend=c("Belief certainty", "Multiple-choice"),
       col=c("black", "gray45"), lty=c(1,1), lwd=4, cex=2, bty="n") 

# TIF
areainfo(nomMC2, c(-3,3)) # 7.23
areainfo(nomTF2, c(-3,3)) # 11.35

# # -----------------------------------------------------------------------
### TIF STUDY 3 
# # -----------------------------------------------------------------------
# The TIF is the sum of the item information functions (IIF) in a test, 
# and indicates the precision of measurement that can be achieved with the 
# test at any value of the latent variable
Theta3 <- matrix(seq(-3,3,.05))
plot(Theta3, testinfo(nomMC3, Theta3), type = 'l', main = 'Study 3 (Bovitz)', 
     ylab = "", xlab=substitute(paste('θ')), ylim=c(0,12), lwd=3, col = 'gray45', lty=1,
     cex.main=2, cex.axis=2, cex.lab=3, frame.plot=F, mgp = c(3.5, 0.5, -0.5))
lines(Theta3, testinfo(nomTF3, Theta3), col = 'black', lty=1, lwd=3)
mtext(side=2, line=2.5, "Information", cex=1.5, font=1)
legend("topright", legend=c("Belief certainty", "Multiple-choice"),
       col=c("black", "gray45"), lty=c(1,1), lwd=4, cex=2, bty="n")

# Area under the curves.
areainfo(nomMC3, c(-3,3)) # 3.49
areainfo(nomTF3, c(-3,3)) # 9.81

dev.off()

# END
print("Code logged successfully")
